nz.map <- ne_countries(country = 'New Zealand', returnclass = "sf", scale = 'medium')
nz.crs <- crs(nz.map)
nz_cropped <- st_crop(nz.map, xmin = 130, xmax = 180,
ymin = -70, ymax = -33)
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
ggplot() + geom_sf(data = nz_cropped) + theme_bw()
# Get the coastline from the rnaturalearth package for comparison
nz <- ne_countries(scale = "medium", returnclass = "sf", country = "New Zealand") %>%
st_transform(crs = proj_nzsf()) %>%
st_crop(get_statistical_areas(area = "EEZ"))
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
Let’s look at the raw catalogue
head ../../../data/NewZealand/hypnzw23R2001_2011.xyze
## YR MODY HRMN SEC LATITUDE LONGITUDE DEPTH MAG NO RMSRES x y GAP DMIN RZDM NP NS SE_OT SE_H SE_Z Q
## 2001 0101 0048 35.29 -39.5230 175.6982 -0.42 2.00 6 0.22 26.16 -342.11 351 33 0.0 3 3 0.24 16.17 23.37 E
## 2001 0101 0153 46.26 -42.2526 172.6404 7.57 2.60 10 0.07 18.33 55.54 172 59 0.1 5 5 0.06 1.00 2.00 C
## 2001 0101 0224 43.34 -38.0645 178.6776 35.40 2.80 6 0.02 133.69 -628.27 290 38 0.9 3 3 0.07 2.00 4.00 E
## 2001 0101 0301 37.34 -38.0325 178.7334 29.31 3.10 7 0.05 135.52 -634.06 255 42 0.7 4 3 0.14 1.50 3.00 D
## 2001 0101 0311 17.49 -38.0626 178.7789 30.95 2.90 7 0.07 140.73 -633.74 265 46 0.7 4 3 0.10 1.50 3.00 D
## 2001 0101 0318 23.15 -41.2742 172.6529 151.57 3.50 22 0.11 -51.04 -28.06 112 53 2.9 16 6 0.08 0.60 0.85 A
## 2001 0101 0403 21.72 -39.2606 175.4920 6.87 1.30 8 0.11 -5.86 -353.93 302 5 1.4 4 4 0.12 2.00 4.00 E
## 2001 0101 0437 9.63 -40.2849 173.2374 187.26 3.20 4 0.52 -83.92 -143.86 253 63 3.0 2 2 0.10 10.79 20.00 E
## 2001 0101 0637 19.56 -39.2356 175.5181 11.32 2.00 8 0.09 -5.82 -357.51 186 4 2.8 7 1 0.16 1.00 2.00 C
To parse the datetimes, I first pre-process the catalogue to split the MODY and HRMN columns. This makes it easier to parse later on dealine with these as integers etc. Alternatively, I could have loaded as strings.
Now we can load the catalogue in R and parse the date-times
# Load the reformatted catalogue
df_cat <- read.table("../../../data/NewZealand/hypnzw23R2001_2011_MNreformat.xyze", header=TRUE, sep="", col.names=c("YR","MO","DY","HR","MN","SEC","LATITUDE","LONGITUDE","DEPTH","MAG","NO","RMSRES","x","y","GAP","DMIN","RZDM","NP","NS","SE_OT","SE_H","SE_Z","Q"), fill=TRUE)
# add event number
df_cat$event_num <- seq.int(nrow(df_cat))
# parse date-times and set appropriate digital second resoluton
op <- options(digits.secs=2)
options(op)
df_cat$date <- ymd( paste(df_cat$YR, df_cat$MO, df_cat$DY, sep="-") )
df_cat$dateTime <- ymd_hms( paste(df_cat$YR,"-", df_cat$MO,"-", df_cat$DY," " , df_cat$HR,"-",df_cat$MN,"-" ,df_cat$SEC, sep="") )
## Warning: 11 failed to parse.
# sp version of the catalogue
df_cat.sp <- data.frame(df_cat)
sp::coordinates(df_cat.sp) <- c("LONGITUDE", "LATITUDE")
crs_wgs84 <- CRS(SRS_string='EPSG:4326')
proj4string(df_cat.sp) <- crs_wgs84
# sf version of the catalogue
df_cat.sf<- st_as_sf(df_cat.sp)
head(df_cat)
## YR MO DY HR MN SEC LATITUDE LONGITUDE DEPTH MAG NO RMSRES x y
## 1 2001 1 1 0 48 35.29 -39.5230 175.6982 -0.42 2.0 6 0.22 26.16 -342.11
## 2 2001 1 1 1 53 46.26 -42.2526 172.6404 7.57 2.6 10 0.07 18.33 55.54
## 3 2001 1 1 2 24 43.34 -38.0645 178.6776 35.40 2.8 6 0.02 133.69 -628.27
## 4 2001 1 1 3 1 37.34 -38.0325 178.7334 29.31 3.1 7 0.05 135.52 -634.06
## 5 2001 1 1 3 11 17.49 -38.0626 178.7789 30.95 2.9 7 0.07 140.73 -633.74
## 6 2001 1 1 3 18 23.15 -41.2742 172.6529 151.57 3.5 22 0.11 -51.04 -28.06
## GAP DMIN RZDM NP NS SE_OT SE_H SE_Z Q event_num date
## 1 351 33 0.0 3 3 0.24 16.17 23.37 E 1 2001-01-01
## 2 172 59 0.1 5 5 0.06 1.00 2.00 C 2 2001-01-01
## 3 290 38 0.9 3 3 0.07 2.00 4.00 E 3 2001-01-01
## 4 255 42 0.7 4 3 0.14 1.50 3.00 D 4 2001-01-01
## 5 265 46 0.7 4 3 0.10 1.50 3.00 D 5 2001-01-01
## 6 112 53 2.9 16 6 0.08 0.60 0.85 A 6 2001-01-01
## dateTime
## 1 2001-01-01 00:48:35
## 2 2001-01-01 01:53:46
## 3 2001-01-01 02:24:43
## 4 2001-01-01 03:01:37
## 5 2001-01-01 03:11:17
## 6 2001-01-01 03:18:23
ggplot(df_cat, aes( x=DEPTH )) +
geom_histogram( binwidth=1 ) +
coord_flip (xlim = c(350, 0) ) +
ggtitle("Hight res. histogram of event depths for depth inversion aretfacts")
## Plot of event magnitudes using natural time
ggplot(df_cat, aes(x=event_num, y=MAG)) + geom_point(size = 0.1)
ggplot(df_cat, aes(x=dateTime, y=MAG)) +
geom_point(size = 0.1) +
ggtitle("New Zealand magnitude time series")
## Warning: Removed 11 rows containing missing values (geom_point).
# Filtered for M>4
ggplot(df_cat[df_cat$MAG>4,], aes(x=dateTime, y=MAG)) +
geom_point(size = 0.1) +
ggtitle("New Zealand magnitude timeseries for M>4")
minMag <- 1
maxMag <- max(df_cat$MAG)
mags <- df_cat[df_cat$MAG>minMag,]$MAG
tmp <- hist(mags, breaks=seq(minMag-0.05,maxMag+0.1,0.1), plot=FALSE)
N.counts <- length( tmp$counts)
tmp$cumulativeCounts <- cumsum(tmp$counts[N.counts:1])[N.counts:1]
m.min <- 4
bin_m.min <- which(tmp$mids==m.min)
freq_m.min <- tmp$counts[bin_m.min]
b <- 1.1
x <- tmp$mids
y <- freq_m.min * 10^(-b*(x-m.min))
y.cum <- tmp$cumulativeCounts[bin_m.min] * 10^(-b*(x-m.min))
ggplot() +
geom_point( aes(x=tmp$mids, y=tmp$counts) ) +
geom_point( aes(x=tmp$mids, y=tmp$cumulativeCounts) , color='red', pch="+", size=2) +
scale_y_log10() +
ggtitle(paste("Frequency-magnitude plot with arbitary GR dist: b =", b)) +
xlab("Magnitude") +
ylab("log10(Frequency)") +
geom_line(aes(x=x, y=y)) +
geom_line(aes(x=x, y=y.cum), color='red') +
geom_vline( xintercept=m.min, lty=2 )
## Warning: Transformation introduced infinite values in continuous y-axis
b.stability.list <- c()
b.error.list <- c()
m.mean <- c()
max.index.x <- length(x)-5
for( i in 1:max.index.x ){
mag.threshold <- x[i]
N.events <- length(mags[mags > mag.threshold])
m.mean <- mean( mags[mags > mag.threshold] )
m.sd <- sd(mags[mags > mag.threshold])
b.stability.list[i] <- log10( exp(1) ) / (m.mean - mag.threshold - 0.05)
b.error.list[i] <- 2.3 * b.stability.list[i]^2 * m.sd / (sqrt(N.events*(N.events-1)))
}
ggplot() +
geom_line( aes(x=x[1:max.index.x], y=b.stability.list) ) +
geom_line( aes(x=x[1:max.index.x], y=b.stability.list+b.error.list), color=2, lty=2 ) +
geom_line( aes(x=x[1:max.index.x], y=b.stability.list-b.error.list), color=2, lty=2 ) +
xlab("Magnitude threshold") +
ylab("Aki b-value estimate") +
geom_hline(yintercept = 1, lty=3) +
ggtitle("b-value stability plot for New Zealand")
ggplot() +
geom_hex(data = df_cat[df_cat$MAG>4,], aes(x = LONGITUDE, y = LATITUDE), bins = 50) +
scale_fill_continuous(type = "viridis") +
geom_sf(data = nz_cropped, fill=alpha("lightgrey", 0), color = 'orange', size=0.2) +
ggtitle("Density plot for M>4 events") +
theme_bw()
ggplot() +
geom_sf(data = df_cat.sf[df_cat$MAG>4,], size = 0.05) +
geom_sf(data = nz_cropped, fill=alpha("lightgrey", 0), color = 'green') +
geom_sf(data = df_cat.sf[df_cat$MAG>6,], size = 0.5, color='orange') +
geom_sf(data = df_cat.sf[df_cat$MAG>7,], size = 0.5, color='red') +
ggtitle("Map of event locations")
ggplot() +
geom_sf(data = nz, fill = alpha("lightgrey", 0), colour = "green") +
geom_sf(data = df_cat.sf[df_cat$MAG>4,], size = 0.05) +
geom_sf(data = nz, fill = alpha("lightgrey", 0), colour = "green") +
annotation_north_arrow(location = "tl", which_north = "grid", style = north_arrow_orienteering, height=unit(1., "cm"), width=unit(1., "cm")) +
annotation_scale(location = "br", unit_category = "metric") +
geom_sf(data = df_cat.sf[df_cat$MAG>6,], size = 0.5, color='orange') +
geom_sf(data = df_cat.sf[df_cat$MAG>7,], size = 0.5, color='red', show.legend=points) +
ggtitle("Map of event locations")
ggplot() +
geom_point(data=df_cat[df_cat$MAG>4,], aes(dateTime, LATITUDE), size=0.1) +
geom_point(data=df_cat[df_cat$MAG>6,], aes(dateTime, LATITUDE), size=1.2, color='orange') +
geom_point(data=df_cat[df_cat$MAG>7,], aes(dateTime, LATITUDE), size=1.5, color='red') +
ggtitle("New Zealand latitude-time plot")
#eventDate <- ymd_hms( '2009-07-15 09:22:29' )
latLims <- c( -42, -47)
minMAG <- 4
# Subset the main catalogue
df_cat.south <- df_cat[df_cat$MAG >= minMAG, ]
df_cat.south <- df_cat.south[ (df_cat.south$LATITUDE<latLims[1]), ]
df_cat.south <- df_cat.south[ (df_cat.south$LATITUDE>latLims[2]), ]
#df_cat.subset <- df_cat.subset[ (df_cat.subset$LONGITUDE>longLims[1]), ]
#df_cat.subset <- df_cat.subset[ (df_cat.subset$LONGITUDE<longLims[2]), ]
endDate <- max(df_cat.south$dateTime)
startDate <- min(df_cat.south$dateTime)
longLims <- c( min(df_cat.south$LONGITUDE), max(df_cat.south$LONGITUDE) )
head(df_cat.south)
## YR MO DY HR MN SEC LATITUDE LONGITUDE DEPTH MAG NO RMSRES x
## 163 2001 1 4 20 9 58.04 -45.0396 167.3916 91.06 4.0 18 0.14 -88.59
## 796 2001 1 19 6 54 23.42 -43.6561 169.6835 5.00 4.0 22 0.14 -60.24
## 1067 2001 1 22 20 41 24.25 -44.9499 167.4568 101.89 4.6 21 0.09 -91.82
## 1496 2001 1 31 0 57 22.26 -43.4153 179.7593 36.00 4.0 43 0.47 557.52
## 1552 2001 2 1 5 2 29.32 -44.5924 168.1326 76.88 4.0 21 0.13 -80.52
## 2190 2001 2 12 3 6 43.91 -46.7453 165.1765 12.75 4.3 7 0.35 -76.92
## y GAP DMIN RZDM NP NS SE_OT SE_H SE_Z Q event_num date
## 163 569.18 187 52 1.8 12 6 0.08 1.00 2.00 C 163 2001-01-04
## 796 331.81 139 35 0.1 17 5 0.02 1.00 5.00 C 796 2001-01-19
## 1067 558.42 191 65 1.6 13 8 0.06 1.00 2.00 C 1067 2001-01-22
## 1496 -199.50 305 433 0.1 39 4 0.32 2.38 4.00 E 1496 2001-01-31
## 1552 492.64 209 86 0.9 14 7 0.05 1.00 2.00 C 1552 2001-02-01
## 2190 825.47 315 244 0.1 6 1 1.03 21.36 26.15 E 2190 2001-02-12
## dateTime
## 163 2001-01-04 20:09:58
## 796 2001-01-19 06:54:23
## 1067 2001-01-22 20:41:24
## 1496 2001-01-31 00:57:22
## 1552 2001-02-01 05:02:29
## 2190 2001-02-12 03:06:43
ggplot() +
geom_sf(data = df_cat.sf[df_cat$MAG>4,], size = 0.05) +
geom_sf(data = nz_cropped, fill=alpha("lightgrey", 0), color = 'green') +
geom_sf(data = df_cat.sf[df_cat$MAG>6,], size = 0.5, color='orange') +
geom_sf(data = df_cat.sf[df_cat$MAG>7,], size = 0.5, color='red') +
ggtitle("Map of event locations")+
geom_rect(aes(xmin=longLims[1], xmax=longLims[2], ymin=latLims[1], ymax=latLims[2]), color="blue",fill=NA)
ggplot() +
geom_point(data=df_cat[df_cat$MAG>4,], aes(dateTime, LATITUDE), size=0.1) +
geom_point(data=df_cat[df_cat$MAG>6,], aes(dateTime, LATITUDE), size=1.2, color='orange') +
geom_point(data=df_cat[df_cat$MAG>7,], aes(dateTime, LATITUDE), size=1.5, color='red') +
ggtitle("New Zealand latitude-time plot with Southern Island highlighted") +
geom_rect(aes(xmin = as.POSIXct(ymd("2001-1-1")), xmax = as.POSIXct(ymd("2012-1-1")), ymin = latLims[1], ymax = latLims[2]), alpha = 0.4, fill='blue')
ggplot(df_cat.south, aes(x=dateTime, y=MAG)) +
geom_point(size = 0.1) +
ggtitle("Southern Island magnitude timeseries for M>4") +
geom_rect( aes(xmin = as.POSIXct(startDate), xmax = as.POSIXct(endDate), ymin = minMAG, ymax = max(df_cat.south$MAG+0.2)), alpha = 0.4, fill=NA, color="blue" )
## Warning: Use of `df_cat.south$MAG` is discouraged. Use `MAG` instead.
minMag <- 4
maxMag <- max(df_cat.south$MAG)
mags <- df_cat.south[df_cat.south$MAG>=minMag,]$MAG
tmp <- hist(mags, breaks=seq(minMag-0.05,maxMag+0.1,0.1), plot=FALSE)
N.counts <- length( tmp$counts)
tmp$cumulativeCounts <- cumsum(tmp$counts[N.counts:1])[N.counts:1]
m.min <- 4
bin_m.min <- which(tmp$mids==m.min)
freq_m.min <- tmp$counts[bin_m.min]
b <- 1.1
x <- tmp$mids
y <- freq_m.min * 10^(-b*(x-m.min))
y.cum <- tmp$cumulativeCounts[bin_m.min] * 10^(-b*(x-m.min))
ggplot() +
geom_point( aes(x=tmp$mids, y=tmp$counts) ) +
geom_point( aes(x=tmp$mids, y=tmp$cumulativeCounts) , color='red', pch="+") +
scale_y_log10() +
ggtitle(paste("Frequency-magnitude plot with arbitary GR dist: b =", b)) +
xlab("Magnitude") +
ylab("log10(Frequency)") +
geom_line(aes(x=x, y=y)) +
geom_line(aes(x=x, y=y.cum), color='red') +
geom_vline( xintercept=m.min, lty=2 )
## Warning: Transformation introduced infinite values in continuous y-axis
b.stability.list <- c()
b.error.list <- c()
m.mean <- c()
max.index.x <- length(x)-5
for( i in 1:max.index.x ){
mag.threshold <- x[i]
m.mean[i] <- mean( mags[mags > mag.threshold], na.rm=TRUE )
b.stability.list[i] <- log10( exp(1) ) / (m.mean[i] - mag.threshold - 0.05)
b.error.list[i] <- 2.3 * b.stability.list[i]^2 * sd(mags[mags > mag.threshold]) / (sqrt(N.events*(N.events-1)))
}
ggplot() +
geom_line( aes(x=x[1:max.index.x], y=b.stability.list) ) +
geom_line( aes(x=x[1:max.index.x], y=b.stability.list+b.error.list), color=2, lty=2 ) +
geom_line( aes(x=x[1:max.index.x], y=b.stability.list-b.error.list), color=2, lty=2 ) +
xlab("Magnitude threshold") +
ylab("Aki b-value estimate") +
geom_hline(yintercept = 1, lty=3) +
ggtitle("b-value stability plot for Southern Island")
#eventDate <- ymd_hms( '2009-07-15 09:22:29' )
latLims <- c( -34, -41)
minMAG <- 4
# Subset the main catalogue
df_cat.north <- df_cat[df_cat$MAG >= minMAG, ]
df_cat.north <- df_cat.north[ (df_cat.north$LATITUDE>latLims[2]), ]
#df_cat.north <- df_cat.north[ (df_cat.north$LATITUDE<latLims[2]), ]
#df_cat.subset <- df_cat.subset[ (df_cat.subset$LONGITUDE>longLims[1]), ]
#df_cat.subset <- df_cat.subset[ (df_cat.subset$LONGITUDE<longLims[2]), ]
endDate <- max(df_cat.north$dateTime)
startDate <- min(df_cat.north$dateTime)
longLims <- c( min(df_cat.north$LONGITUDE), max(df_cat.north$LONGITUDE) )
head(df_cat.north)
## YR MO DY HR MN SEC LATITUDE LONGITUDE DEPTH MAG NO RMSRES x
## 14 2001 1 1 11 58 50.85 -38.0754 178.7195 27.32 4.0 25 0.08 137.38
## 39 2001 1 2 5 34 42.69 -39.0541 176.8003 38.82 4.5 46 0.29 69.16
## 42 2001 1 2 6 22 16.54 -40.5215 175.6788 37.57 4.4 51 0.14 93.54
## 150 2001 1 4 15 46 53.47 -38.0652 178.7116 33.04 4.8 46 0.15 136.14
## 180 2001 1 5 5 30 45.62 -38.4291 175.6323 179.26 4.1 31 0.20 -53.64
## 203 2001 1 5 17 4 38.86 -38.6254 175.1182 234.35 4.1 39 0.15 -75.15
## y GAP DMIN RZDM NP NS SE_OT SE_H SE_Z Q event_num date
## 14 -629.49 252 31 0.9 18 7 0.06 1.50 3.00 D 14 2001-01-01
## 39 -441.69 82 49 0.8 41 5 0.08 0.41 1.79 B 39 2001-01-02
## 42 -254.02 74 18 2.1 40 11 0.02 0.24 0.32 A 42 2001-01-02
## 150 -629.99 249 40 0.8 36 10 0.08 1.50 3.00 D 150 2001-01-04
## 180 -433.87 186 64 2.8 24 7 0.12 1.00 2.00 C 180 2001-01-05
## 203 -388.89 144 35 6.7 31 8 0.12 0.66 1.03 B 203 2001-01-05
## dateTime
## 14 2001-01-01 11:58:50
## 39 2001-01-02 05:34:42
## 42 2001-01-02 06:22:16
## 150 2001-01-04 15:46:53
## 180 2001-01-05 05:30:45
## 203 2001-01-05 17:04:38
ggplot() +
geom_sf(data = df_cat.sf[df_cat$MAG>4,], size = 0.05) +
geom_sf(data = nz_cropped, fill=alpha("lightgrey", 0), color = 'green') +
geom_sf(data = df_cat.sf[df_cat$MAG>6,], size = 0.5, color='orange') +
geom_sf(data = df_cat.sf[df_cat$MAG>7,], size = 0.5, color='red') +
ggtitle("Map of event locations")+
geom_rect(aes(xmin=longLims[1], xmax=longLims[2], ymin=latLims[1], ymax=latLims[2]), color="blue",fill=NA)
ggplot() +
geom_point(data=df_cat[df_cat$MAG>4,], aes(dateTime, LATITUDE), size=0.1) +
geom_point(data=df_cat[df_cat$MAG>6,], aes(dateTime, LATITUDE), size=1.2, color='orange') +
geom_point(data=df_cat[df_cat$MAG>7,], aes(dateTime, LATITUDE), size=1.5, color='red') +
ggtitle("New Zealand latitude-time plot with Northern Island highlighted") +
geom_rect(aes(xmin = as.POSIXct(ymd("2001-1-1")), xmax = as.POSIXct(ymd("2012-1-1")), ymin = latLims[1], ymax = latLims[2]), alpha = 0.4, fill='blue')
ggplot() +
geom_point(data=df_cat.north[df_cat.north$MAG>4,], aes(dateTime, LATITUDE), size=0.1) +
geom_point(data=df_cat.north[df_cat.north$MAG>6,], aes(dateTime, LATITUDE), size=1.2, color='orange') +
geom_point(data=df_cat.north[df_cat.north$MAG>7,], aes(dateTime, LATITUDE), size=1.5, color='red') +
ggtitle("New Zealand latitude-time plot with Northern Island highlighted") +
geom_rect(aes(xmin = as.POSIXct(ymd("2001-1-1")), xmax = as.POSIXct(ymd("2012-1-1")), ymin = latLims[1], ymax = latLims[2]), alpha = 0.4, fill='blue')
ggplot(df_cat.north, aes(x=dateTime, y=MAG)) +
geom_point(size = 0.1) +
ggtitle("Northern Island magnitude timeseries for M>4") +
geom_rect( aes(xmin = as.POSIXct(startDate), xmax = as.POSIXct(endDate), ymin = minMAG, ymax = max(df_cat.north$MAG+0.2)), alpha = 0.4, fill=NA, color="blue" )
## Warning: Use of `df_cat.north$MAG` is discouraged. Use `MAG` instead.
minMag <- 4
maxMag <- max(df_cat.north$MAG)
mags <- df_cat.north[df_cat.north$MAG>=minMag,]$MAG
tmp <- hist(mags, breaks=seq(minMag-0.05,maxMag+0.1,0.1), plot=FALSE)
N.counts <- length( tmp$counts)
tmp$cumulativeCounts <- cumsum(tmp$counts[N.counts:1])[N.counts:1]
m.min <- 4
bin_m.min <- which(tmp$mids==m.min)
freq_m.min <- tmp$counts[bin_m.min]
b <- 1.1
x <- tmp$mids
y <- freq_m.min * 10^(-b*(x-m.min))
y.cum <- tmp$cumulativeCounts[bin_m.min] * 10^(-b*(x-m.min))
ggplot() +
geom_point( aes(x=tmp$mids, y=tmp$counts) ) +
geom_point( aes(x=tmp$mids, y=tmp$cumulativeCounts) , color='red', pch="+") +
scale_y_log10() +
ggtitle(paste("Frequency-magnitude plot with arbitary GR dist: b =", b)) +
xlab("Magnitude") +
ylab("log10(Frequency)") +
geom_line(aes(x=x, y=y)) +
geom_line(aes(x=x, y=y.cum), color='red') +
geom_vline( xintercept=m.min, lty=2 )
## Warning: Transformation introduced infinite values in continuous y-axis
b.stability.list <- c()
b.error.list <- c()
m.mean <- c()
max.index.x <- length(x)-5
for( i in 1:max.index.x ){
mag.threshold <- x[i]
m.mean[i] <- mean( mags[mags > mag.threshold], na.rm=TRUE )
b.stability.list[i] <- log10( exp(1) ) / (m.mean[i] - mag.threshold - 0.05)
b.error.list[i] <- 2.3 * b.stability.list[i]^2 * sd(mags[mags > mag.threshold]) / (sqrt(N.events*(N.events-1)))
}
ggplot() +
geom_line( aes(x=x[1:max.index.x], y=b.stability.list) ) +
geom_line( aes(x=x[1:max.index.x], y=b.stability.list+b.error.list), color=2, lty=2 ) +
geom_line( aes(x=x[1:max.index.x], y=b.stability.list-b.error.list), color=2, lty=2 ) +
xlab("Magnitude threshold") +
ylab("Aki b-value estimate") +
geom_hline(yintercept = 1, lty=3) +
ggtitle("b-value stability plot for Northern Island")
UTC time 2009-07-15 09:22:29 ISC event 15157724 USGS-ANSS ComCat Local date 15 July 2009 Local time 9:22 pm (NZST) Magnitude 7.8 Mw Depth 12.0 km (7.5 mi) Epicentre 45.762°S 166.562°ECoordinates: 45.762°S 166.562°E Areas affected New Zealand Max. intensity VI (Strong) Tsunami 2.3 m (7 ft 7 in) Aftershocks >100 Casualties 0
eventDate <- ymd_hms( '2009-07-15 09:22:29' )
endDate <- eventDate + days(100)
startDate <- eventDate - days(50)
deltaLat <- 2.4
latLims <- c( -45.762-deltaLat, -45.762+deltaLat)
longLims <- c( 166.562-deltaLat, 166.562+deltaLat)
minMAG <- 4
# Subset the main catalogue
df_cat.subset <- df_cat[df_cat$MAG >= minMAG, ]
df_cat.subset <- df_cat.subset[ (df_cat.subset$LATITUDE>latLims[1]), ]
df_cat.subset <- df_cat.subset[ (df_cat.subset$LATITUDE<latLims[2]), ]
df_cat.subset <- df_cat.subset[ (df_cat.subset$LONGITUDE>longLims[1]), ]
df_cat.subset <- df_cat.subset[ (df_cat.subset$LONGITUDE<longLims[2]), ]
head(df_cat.subset)
## YR MO DY HR MN SEC LATITUDE LONGITUDE DEPTH MAG NO RMSRES x
## 163 2001 1 4 20 9 58.04 -45.0396 167.3916 91.06 4.0 18 0.14 -88.59
## 1067 2001 1 22 20 41 24.25 -44.9499 167.4568 101.89 4.6 21 0.09 -91.82
## 1552 2001 2 1 5 2 29.32 -44.5924 168.1326 76.88 4.0 21 0.13 -80.52
## 2190 2001 2 12 3 6 43.91 -46.7453 165.1765 12.75 4.3 7 0.35 -76.92
## 2194 2001 2 12 4 34 17.26 -46.8243 165.3218 30.00 5.0 11 0.38 -62.82
## 3122 2001 2 25 20 1 28.97 -44.6792 168.1733 11.52 4.3 23 0.14 -71.54
## y GAP DMIN RZDM NP NS SE_OT SE_H SE_Z Q event_num date
## 163 569.18 187 52 1.8 12 6 0.08 1.00 2.00 C 163 2001-01-04
## 1067 558.42 191 65 1.6 13 8 0.06 1.00 2.00 C 1067 2001-01-22
## 1552 492.64 209 86 0.9 14 7 0.05 1.00 2.00 C 1552 2001-02-01
## 2190 825.47 315 244 0.1 6 1 1.03 21.36 26.15 E 2190 2001-02-12
## 2194 823.67 319 215 0.1 7 4 0.71 8.41 16.22 E 2194 2001-02-12
## 3122 497.43 129 20 0.6 17 6 0.02 0.25 0.64 A 3122 2001-02-25
## dateTime
## 163 2001-01-04 20:09:58
## 1067 2001-01-22 20:41:24
## 1552 2001-02-01 05:02:29
## 2190 2001-02-12 03:06:43
## 2194 2001-02-12 04:34:17
## 3122 2001-02-25 20:01:28
ggplot() +
geom_sf(data = df_cat.sf[df_cat$MAG>4,], size = 0.05) +
geom_sf(data = nz_cropped, fill=alpha("lightgrey", 0), color = 'green') +
geom_sf(data = df_cat.sf[df_cat$MAG>6,], size = 0.5, color='orange') +
geom_sf(data = df_cat.sf[df_cat$MAG>7,], size = 0.5, color='red') +
ggtitle("Map of event locations")+
geom_rect(aes(xmin=longLims[1], xmax=longLims[2], ymin=latLims[1], ymax=latLims[2]), color="blue",fill=NA)
ggplot() +
geom_point(data=df_cat[df_cat$MAG>4,], aes(dateTime, LATITUDE), size=0.1) +
geom_point(data=df_cat[df_cat$MAG>6,], aes(dateTime, LATITUDE), size=1.2, color='orange') +
geom_point(data=df_cat[df_cat$MAG>7,], aes(dateTime, LATITUDE), size=1.5, color='red') +
ggtitle("New Zealand latitude-time plot") +
geom_rect(aes(xmin = as.POSIXct(ymd("2001-1-1")), xmax = as.POSIXct(ymd("2012-1-1")), ymin = latLims[1], ymax = latLims[2]), alpha = 0.4, fill='blue')
ggplot() +
geom_point(data=df_cat[df_cat$MAG>4,], aes(dateTime, LATITUDE), size=0.1) +
geom_point(data=df_cat[df_cat$MAG>6,], aes(dateTime, LATITUDE), size=1.2, color='orange') +
geom_point(data=df_cat[df_cat$MAG>7,], aes(dateTime, LATITUDE), size=1.5, color='red') +
ggtitle("New Zealand latitude-time plot") +
geom_rect(aes(xmin = as.POSIXct(startDate), xmax = as.POSIXct(endDate), ymin = latLims[1], ymax = latLims[2]), alpha = 0.4, fill='blue')
ggplot(df_cat.subset, aes(x=dateTime, y=MAG)) +
geom_point(size = 0.1) +
ggtitle("New Zealand magnitude timeseries for M>4") +
geom_rect( aes(xmin = as.POSIXct(startDate), xmax = as.POSIXct(endDate), ymin = minMAG, ymax = max(df_cat.subset$MAG+0.2)), alpha = 0.4, fill=NA, color="blue" )
## Warning: Use of `df_cat.subset$MAG` is discouraged. Use `MAG` instead.
minMag <- 4
maxMag <- max(df_cat.subset$MAG)
mags <- df_cat.subset[df_cat.subset$MAG>=minMag,]$MAG
tmp <- hist(mags, breaks=seq(minMag-0.05,maxMag+0.1,0.1), plot=FALSE)
N.counts <- length( tmp$counts)
tmp$cumulativeCounts <- cumsum(tmp$counts[N.counts:1])[N.counts:1]
m.min <- 4
bin_m.min <- which(tmp$mids==m.min)
freq_m.min <- tmp$counts[bin_m.min]
b <- 1.1
x <- tmp$mids
y <- freq_m.min * 10^(-b*(x-m.min))
y.cum <- tmp$cumulativeCounts[bin_m.min] * 10^(-b*(x-m.min))
ggplot() +
geom_point( aes(x=tmp$mids, y=tmp$counts) ) +
geom_point( aes(x=tmp$mids, y=tmp$cumulativeCounts) , color='red', pch="+") +
scale_y_log10() +
ggtitle(paste("Frequency-magnitude plot with arbitary GR dist: b =", b)) +
xlab("Magnitude") +
ylab("log10(Frequency)") +
geom_line(aes(x=x, y=y)) +
geom_line(aes(x=x, y=y.cum), color='red') +
geom_vline( xintercept=m.min, lty=2 )
## Warning: Transformation introduced infinite values in continuous y-axis
b.stability.list <- c()
b.error.list <- c()
m.mean <- c()
max.index.x <- length(x)-5
for( i in 1:max.index.x ){
mag.threshold <- x[i]
m.mean[i] <- mean( mags[mags > mag.threshold], na.rm=TRUE )
b.stability.list[i] <- log10( exp(1) ) / (m.mean[i] - mag.threshold - 0.05)
b.error.list[i] <- 2.3 * b.stability.list[i]^2 * sd(mags[mags > mag.threshold]) / (sqrt(N.events*(N.events-1)))
}
ggplot() +
geom_line( aes(x=x[1:max.index.x], y=b.stability.list) ) +
geom_line( aes(x=x[1:max.index.x], y=b.stability.list+b.error.list), color=2, lty=2 ) +
geom_line( aes(x=x[1:max.index.x], y=b.stability.list-b.error.list), color=2, lty=2 ) +
xlab("Magnitude threshold") +
ylab("Aki b-value estimate") +
geom_hline(yintercept = 1, lty=3) +
ggtitle("b-value stability plot for catalogue subset")